function [] = runSummarizeOutcomeVariable(runId)

%%% Add paths
addpath('../general_funcs');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% Get input / output paths
config = getConfig_jointProcess(runId);

% Get config_consumerAdoptions and config_providerAdoptions: necessary input for make_Xdynamic_xxx functions
indivConfigs.consumerAdoptions  = getConfig_consumerAdoptions(config.runId_consumerAdoptions);
indivConfigs.providerAdoptions  = getConfig_providerAdoptions(config.runId_providerAdoptions);
indivConfigs.providerExits      = getConfig_providerExits(config.runId_providerExits);

indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;

% Load the data
load(sprintf('%s/runData.mat', indivRunFolders.consumerAdoptions));
mydata.consumerAdoptions = data; clear data;
load(sprintf('%s/runData.mat', indivRunFolders.providerAdoptions));
mydata.providerAdoptions = data; clear data;
load(sprintf('%s/runData.mat', indivRunFolders.providerExits));
mydata.providerExits = data; clear data;
data = mydata;

disp('===== Distribution of spell durations =====')
spellDurations = 24*60*data.consumerAdoptions.spellDurations; % in minutes
summary_func(spellDurations)

% Make plots of spell durations (histogram and percentile plot)
makeHistogram_and_percentilePlot(spellDurations, config.runFolder, 'spellDurations', 'Spell duration (minutes)', 100, 400);

disp('===== Distribution of distance between consumer and provider across transactions =====')
T = data.consumerAdoptions.dims.dims1{4};
K1 = data.consumerAdoptions.dims.dims1{5};
K2 = data.consumerAdoptions.dims.dims1{6};

% Count number of transactions completed for each (k1,k2) pair (summed up over time)
Y.consumerAdoptions = reshape(data.consumerAdoptions.Y, [T K1*K2]);
sumY1 = sum(Y.consumerAdoptions,1)'; % (K1*K2) x 1

% Load matrix of distances
load('../../data/data_mat/SpatialWeights/GeoDistance.mat');
%	W(k1,k2) gives distance between centroids of cantons k1 and k2 (in kilometers)
W = reshape(W, [K1*K2 1]); % (K1*K2) x 1

% Construct vector of distances during transactions
NumTransactions = sum(sumY1);
transacDistances = zeros(NumTransactions, 1);

idxes = find(sumY1>0);
sumY1 = sumY1(idxes); % NumUnqPairs x 1 (unique pairs with transactions)
W = W(idxes); % NumUnqPairs x 1 (unique pairs with transactions)
offset = 0;
for ii = 1:length(sumY1)
	pairDistance = W(ii);
	pairY = sumY1(ii);
	transacDistances(offset+1:offset+pairY) = pairDistance;
	
	offset = offset + pairY;
end
assert(offset == NumTransactions); % simple check

summary_func(transacDistances)

disp('Percentage of transactions where consumer-provider distance is less than 40km:')
mean(transacDistances < 40)

% Make plots of spell durations (histogram and percentile plot)
makeHistogram_and_percentilePlot(transacDistances, config.runFolder, 'transacDistances', 'Distance consumer-provider (kilometer)', 30, 1200);


disp('===== Distribution of count of consumer adoptions / provider adoptions / provider exits per day and per week =====')
NumSpells = data.consumerAdoptions.dims.dims1{4};
K1 = data.consumerAdoptions.dims.dims1{5};
K2 = data.consumerAdoptions.dims.dims1{6};
K = data.providerAdoptions.dims.dims1{3};

Y.consumerAdoptions = reshape(data.consumerAdoptions.Y, [T K1*K2]); % NumSpells x (K1*K2)
Y.providerAdoptions  = reshape(data.providerAdoptions.Y, [T K]); % NumSpells x K
Y.providerExits      = reshape(data.providerExits.Y, [T K]);     % NumSpells x K

spell2period = data.consumerAdoptions.spell2period; % NumSpells x 1 (values between 1 and NumDays)
NumDays = max(spell2period);

% Aggregate to the daily level (take sum across spells of that day tt)
M = sparse(NumDays, NumSpells);
ind = sub2ind([NumDays NumSpells], spell2period, [1:NumSpells]');
M(ind) = 1;
Ydaily.consumerAdoptions = M * Y.consumerAdoptions; % NumDays x (K1*K2)
Ydaily.providerAdoptions  = M * Y.providerAdoptions; % NumDays x K
Ydaily.providerExits  = M * Y.providerExits; % NumDays x K
clear M ind;

% Aggregate to the weekly level
day2week = ceil([1:double(NumDays)]./7)';
NumWeeks = max(day2week);
M = zeros(NumWeeks, NumDays);
ind = sub2ind([NumWeeks NumDays], day2week, [1:NumDays]');
M(ind) = 1;
Yweekly.consumerAdoptions = M * Ydaily.consumerAdoptions; % NumWeeks x (K1*K2)
Yweekly.providerAdoptions  = M * Ydaily.providerAdoptions; % NumWeeks x K
Yweekly.providerExits  = M * Ydaily.providerExits; % NumWeeks x K
clear M ind;

% Expand as vectors
Ydaily.consumerAdoptions = Ydaily.consumerAdoptions(:); % (NumDays*K1*K2) x 1
Ydaily.providerAdoptions  = Ydaily.providerAdoptions(:);  % (NumDays*K) x 1
Ydaily.providerExits      = Ydaily.providerExits(:);      % (NumDays*K) x 1

Yweekly.consumerAdoptions = Yweekly.consumerAdoptions(:); % (NumWeeks*K1*K2) x 1
Yweekly.providerAdoptions  = Yweekly.providerAdoptions(:);  % (NumWeeks*K) x 1
Yweekly.providerExits      = Yweekly.providerExits(:);      % (NumWeeks*K) x 1


%%% Put together distribution level at daily level
% Look at distribution of consumerAdoptions at daily level
idx = find(Ydaily.consumerAdoptions);
NumZeros = length(Ydaily.consumerAdoptions) - length(idx);
tmp = tabulate(Ydaily.consumerAdoptions(idx));
unqValues = tmp(:,1);
freqs     = tmp(:,2);
unqValues = [0; unqValues];
freqs     = [NumZeros; freqs];
rel_freqs = 100*freqs/sum(freqs);
Mday_consumerAdoptions = [unqValues freqs rel_freqs];
clear idx NumZeros tmp unqValues freqs rel_freqs

% Look at distribution of providerAdoptions at daily level
Mday_providerAdoptions = tabulate(full(Ydaily.providerAdoptions));

% Look at distribution of providerExits at daily level
Mday_providerExits = tabulate(full(Ydaily.providerExits));


%%%% Put together these values
maxValue = max([Mday_consumerAdoptions(:,1); Mday_providerAdoptions(:,1); Mday_providerExits(:,1)]);
myvalues = [0:maxValue]';
NumValues = length(myvalues);
freqs = zeros(NumValues,3);
for vv = 1:NumValues
	val = myvalues(vv);
	% consumerAdoptions (1st column of freqs)
	idx = find(Mday_consumerAdoptions(:,1) == val);
	if ~isempty(idx)
		freqs(vv,1) = Mday_consumerAdoptions(idx,2);
	end
	clear idx;
	
	% providerAdoptions (2nd column of freqs)
	idx = find(Mday_providerAdoptions(:,1) == val);
	if ~isempty(idx)
		freqs(vv,2) = Mday_providerAdoptions(idx,2);
	end
	clear idx;
	
	% providerExits (3rd column of freqs)
	idx = find(Mday_providerExits(:,1) == val);
	if ~isempty(idx)
		freqs(vv,3) = Mday_providerExits(idx,2);
	end
	clear idx;
end
rel_freqs = freqs ./sum(freqs,1);
% Add total
myvaluesStrings = sprintfc('%d', myvalues);
myvaluesStrings = [myvaluesStrings; {'Total'}];
freqs           = [freqs; sum(freqs)];
rel_freqs       = [rel_freqs; sum(rel_freqs)];
% Combine everything into a table and export it to CSV
dailyFreqTable = [cell2table(myvaluesStrings), array2table([freqs rel_freqs])];
dailyFreqTable.Properties.VariableNames = {'Value', 'FreqConsumerAdoptions', 'FreqProviderAdoptions', 'FreqProviderExits', ...
	'RelFreqConsumerAdoptions', 'RelFreqProviderAdoptions', 'RelFreqProviderExits'};
disp(dailyFreqTable);
writetable(dailyFreqTable, sprintf('%s/freqTable_daily.csv', config.runFolder));

clear maxValue myvalues NumValues freqs rel_freqs;



%%% Put together distribution level at weekly level
% Look at distribution at weekly level for consumerAdoptions (complicated because of data size do part of the process myself)
idx = find(Yweekly.consumerAdoptions);
NumZeros = length(Yweekly.consumerAdoptions) - length(idx);
tmp = tabulate(Yweekly.consumerAdoptions(idx));
unqValues = tmp(:,1);
freqs     = tmp(:,2);
unqValues = [0; unqValues];
freqs     = [NumZeros; freqs];
rel_freqs = 100*freqs/sum(freqs);
Mweek_consumerAdoptions = [unqValues freqs rel_freqs];
clear idx NumZeros tmp unqValues freqs rel_freqs;

% Look at distribution of providerAdoptions at weekly level
Mweek_providerAdoptions = tabulate(full(Yweekly.providerAdoptions));

% Look at distribution of providerExits at weekly level
Mweek_providerExits = tabulate(full(Yweekly.providerExits));


%%%% Put together these values
maxValue = max([Mweek_consumerAdoptions(:,1); Mweek_providerAdoptions(:,1); Mweek_providerExits(:,1)]);
myvalues = [0:maxValue]';
NumValues = length(myvalues);
freqs = zeros(NumValues,3);
for vv = 1:NumValues
	val = myvalues(vv);
	% consumerAdoptions (1st column of freqs)
	idx = find(Mweek_consumerAdoptions(:,1) == val);
	if ~isempty(idx)
		freqs(vv,1) = Mweek_consumerAdoptions(idx,2);
	end
	clear idx;
	
	% providerAdoptions (2nd column of freqs)
	idx = find(Mweek_providerAdoptions(:,1) == val);
	if ~isempty(idx)
		freqs(vv,2) = Mweek_providerAdoptions(idx,2);
	end
	clear idx;
	
	% providerExits (3rd column of freqs)
	idx = find(Mweek_providerExits(:,1) == val);
	if ~isempty(idx)
		freqs(vv,3) = Mweek_providerExits(idx,2);
	end
	clear idx;
end
rel_freqs = freqs ./sum(freqs,1);
% Add total
myvaluesStrings = sprintfc('%d', myvalues);
myvaluesStrings = [myvaluesStrings; {'Total'}];
freqs           = [freqs; sum(freqs)];
rel_freqs       = [rel_freqs; sum(rel_freqs)];
% Combine everything into a table and export it to CSV
weeklyFreqTable = [cell2table(myvaluesStrings), array2table([freqs rel_freqs])];
weeklyFreqTable.Properties.VariableNames = {'Value', 'FreqConsumerAdoptions', 'FreqProviderAdoptions', 'FreqProviderExits', ...
	'RelFreqConsumerAdoptions', 'RelFreqProviderAdoptions', 'RelFreqProviderExits'};
disp(weeklyFreqTable);
writetable(weeklyFreqTable, sprintf('%s/freqTable_weekly.csv', config.runFolder));

clear maxValue myvalues NumValues freqs rel_freqs;


end
